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Abstract 

Strong stability preserving (SSP) methods are designed primarily for time integration of 
nonlinear h5rperbolic PDFs, for which fhe permissible SSP sfep size varies from one sfep 
fo fhe nexf. We develop fhe firsf SSP linear mulfisfep mefhods (of order fwo and fhree) 
wifh variable sfep size, and prove fheir opfimalify sfabilify, and convergence. The choice 
of sfep size for mulfisfep SSP mefhods is an inferesfing problem because fhe allowable sfep 
size depends on fhe SSP coefficienf, which in furn depends on fhe chosen sfep sizes. The 
descripfion of fhe mefhods includes an opfimal sfep-size sfrafegy. We prove sharp upper 
bounds on fhe allowable sfep size for explicif SSP linear mulfisfep mefhods and show fhe 
exisfence of mefhods wifh arbifrarily high order of accuracy. The effectiveness of fhe mefhods 
is demonsfrafed fhrough numerical examples. 


1 Introduction 

Strong stability preserving (SSP) linear multistep methods (LMMs) with uniform step size have 
been studied by several authors [Len89; Len91; HRS03; RH05; HR05; Ket09]. In this work, we 
develop the first variable step-size SSP multistep methods. 

The principal area of application of strong stability preserving methods is the integration of 
nonlinear systems of hyperbolic conservation laws. In such applications, the allowable step size 
hmax is usually determined by a CFL-like condition, and in particular is inversely proportional 
to the fastest wave speed. This wave speed may vary significantly during the course of the 
integration, and its variation cannot generally be predicted in advance. Thus a fixed-step-size 
code may be inefficient (if h^ax increases) or may fail completely (if //max decreases). 
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SSP Runge-Kutta methods are used much more widely than SSP LMMs. Indeed, it is difficult 
to find examples of SSP LMMs used in realistic applications; this may be due to the lack of 
a variable step-size (VSS) formulation. Tradeoffs between Runge-Kutta and linear multistep 
methods have been discussed at length elsewhere, but one reason for preferring LMMs over 
their Runge-Kutta counterparts in the context of strong stability preservation stems from the 
recent development of a high-order positivity preserving limiter [ZSIO]. Use of Runge-Kutta 
methods in conjunction with the limiter can lead to order reduction, so linear multistep methods 
are recommended [GKSll; ZSIO]. 

There exist two approaches to variable step-size multistep methods [HNW93]. In the first, 
polynomial interpolation is used to find values at equally spaced points, and then the fixed-step- 
size method is used. In the second, the method coefficients are varied to maintain high order 
accuracy based on the solution values at the given (non-uniform) step intervals. Both approaches 
are problematic for SSP methods; the first, because the interpolation step itself may violate the 
SSP property, and the second, because the SSP step size depends on the method coefficients. 
Herein we pursue the second strategy. 

The main contributions of this work are: 

• sharp bounds on the SSP coefficient of variable step-size SSP LMMs (Theorems 2-4); 

• optimal methods of orders two and three (Sections 3.2 and 3.3); 

• existence of methods of arbitrary order (Theorem 5); 

• analysis of the greedy step-size strategy (Section 4); 

• proof of stability and convergence of the optimal methods (Theorems 8-10). 

The rest of the paper is organized as follows. In Section 2.1, we review the theory of SSP 
LMMs while recognizing that the forward Euler permissible step size may change from step to 
step. The main result. Theorem 1, is a slight refinement of the standard one. In Section 2.2 we 
show how an optimal SSP multistep formula may be chosen at each step, given the sequence of 
previous steps. In Section 2.3 we provide for convenience a description of the two of the simplest 
and most useful methods in this work. In Sections 3.2 and 3.3 we derive and prove the optimality 
of several second- and third-order methods. In Section 4 we investigate the relation between the 
SSP step size, the method coefficients, and the step-size sequence. We develop step-size strategies 
that ensure the SSP property under mild assumptions on the problem. In Section 5 we prove that 
the methods, with the prescribed step-size strategies, are stable and convergent. In Section 6 we 
demonstrate the efficiency of the methods with some numerical examples. Finally, Sections 8 and 
9 contain the proofs of the more technical theorems and lemmas. 

Two topics that might be pursued in the future based on this work are: 

• variable step-size SSP LMMs of order higher than three; 

• variable step-size versions of SSP methods with multiple steps and multiple stages. 
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2 SSP linear multistep methods 

We consider the numerical solution of the initial value problem 


u'{t) = f{u{t)) u{to) = Uo, (1) 

for t G [fo, fo + ^] by an explicit linear multistep method. If a fixed numerical step size h is used, 
the method takes the form 

k-l 

ii„ = {o^jUn-k+j + h^jf{u„_k+j)) n>k. (2) 

/=o 


Here k is the number of steps and is an approximation to the solution u{nh). 

Now let the step size vary from step to step so that + /i„. In order to achieve the 

same order of accuracy, the coefficients cl, j6 must also vary from step to step: 

k-l 

Un — . (3) 

;=0 


At this point, it is helpful to establish the following terminology We use the term multistep 
formula, or ]ust. formula, to refer to a set of coefficients ay (= ay „), fj (= /3y,,) that may be used at 
step n. We use the term multistep method to refer to a full time-stepping algorithm that includes 
a prescription of how to choose a formula (ay „,/3y „) and the step size h,, at step n. 

For 1 < y < k let 


coi := 


-k+j 


hfi 


> 0 


(4) 


denote the step-size ratios and 


' Oq 0, 

i 

Qy : = ^ coi for 1 < j < k. 
- ! = 1 


Note that the values co and O depend on n, but we often suppress that dependence since we are 
considering a single step.^ It is useful to keep in mind that Qy = y if the step size is fixed. Also 
the simple relation -|- 1 will often be used. 

^Our definition of co differs from the typical approach in the literature on variable step-size multistep methods, 
where only ratios of adjacent step sizes are used. The present definition is more convenient in what follows. 


3 





2.1 Strong stability preservation 

We are interested in initial value problems (1) whose solution satisfies a monotonicity condition 

\\u{t + h)\\ < \\u{t)\\ for//> 0, (6) 

where || • || represents any convex functional (for instance, a norm). We assume that / satisfies 
the (stronger) forward Euler condition 

||m + 12/(m)|| < ||m|| for 0 < Iz </jfe(w). (7) 


The discrete monotonicity condition (7) implies the continuous monotonicity condition (6). 

The primary application of SSP methods is in the time integration of nonlinear hyperbolic 
PDEs. In such applications, is proportional to the CEE number 


V = h 


a{u) 

Ax 


( 8 ) 


where a{u) is the largest wave speed appearing in the problem. This speed depends on u. Eor 
instance, in the case of Burgers' equation 

we have a{u) — maxj \ u\. Eor scalar conservation laws like Burgers' equation, it is possible to 
determine a value of hj;^, based on the initial and boundary data, that is valid for all time. But 
for general systems of conservation laws, a{u) can grow in time and so the minimum value of 
cannot be determined without solving the initial value problem. We will often write just 
/zpE for brevity, but the dependence of /zfe on u should be remembered. 

We will develop linear multistep methods (3) that satisfy the discrete monotonicity property 


U„\\ < max(||l/„_;,||,||M„_,t+l||/---/||Wn-l||)- 


( 10 ) 


The class of methods that satisfy (10) whenever / satisfies (7) are known as strong stability 
preserving methods. The most widely used SSP methods are one-step (Runge-Kutta) meth¬ 
ods. When using an SSP multistep method, an SSP Runge-Kutta method can be used to ensure 
monotonicity of the starting values. In the remainder of this work, we focus on conditions for 
monotonicity of subsequent steps (for any given starting values). 

The following theorem refines a well-known result in the literature, by taking into account 
the dependence of /zfe on u. 

Theorem 1. Suppose that f satisfies the forward Euler condition (7) and that the method (3) has non¬ 
negative coefficients aj^„, > 0. Furthermore, suppose that the time step is chosen so that 

0<h„< min (^ h^{u„^k+j)) (11) 

\py,n / 
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for each n, where the ratio ccj^n/ f>j,n understood as +oo if fj^n — 0- Then the solution of the initial value 

problem (1) given by the LMM (3) satisfies the monotonicity condition (10). 

Remark 1. The step-size restriction (11) is also necessary/or monotonicity in the sense that, for any 
method (3), there exists some f and starting values such that the monotonicity condition (10) will be 
violated if the step size (11) is exceeded. 

Remark 2. Even in the case of the fixed-step-size method (2), the theorem above generalizes results in 
the literature that are based on the assumption of a constant %e- T Is natural to implement a step-size 
strategy that uses simply %E(Mn-i) in (11), but this will not give the correct step size in general. On the 
other hand, since /jfe usually varies slowly from one step to the next, is often non-decreasing, and since 
the restriction (11) is often pessimistic, such a strategy will usually work well. 

Since (and hence ^zfe) varies slowly from one step to the next, it seems convenient to 
separate the factors in the upper bound in (11) and consider the sufficient condition 

0<h„< CnPn (12) 


where the SSP coefficient C„ is 
Cn = 


max {r G R+ : — rfj^n > 0} if > 0/ /;> > 0 for all /; 

0 otherwise. 


(13) 


and 


Pn'- min hYE{un-k+i) {n>k). (14) 

0<;<;k-l ^ 

Note that in general the SSP coefficient varies from step to step, since it depends on the method 
coefficients. 


2.2 Optimal SSP formulae 

For a given order p, number of steps k, and previous step-size sequence //„_!, h„^ 2 , ■ ■ ■, we say that 
a multistep formula is optimal if it gives the largest possible SSP coefficient C„ in (12) and satisfies 
the order conditions. In this section we formulate this optimization problem algebraically. The 
linear multistep formula takes the form 


k-l 

nn — <:+;■ 3 “ A:+/)) • (1^) 

M 
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Here we have omitted the subscript n on the method coefficients to simplify the notation. The 
conditions for formula (15) to be consistent of order p are 


k-l 



(16a) 



1 < m < p. 


(16b) 


Let us change variables by introducing 


Sj ■- Oij - r^j. 


Then the SSP coefficient of the formula is just the largest r such that all 5j, /ly are non-negative 
[Len89; Ket09]. In these variables, the order conditions (16) become 


k-l 


Y.{Si+r^,)=l 


(17a) 





(17b) 


1 < m < p. 


We will refer to a formula by the triplet Given p, k, and a set of step-size ratios tOj, the 

formula with the largest SSP coefficient for the next step can be obtained by finding the largest r 
such that (17) has a non-negative solution ci,/3 > 0. This could be done following the approach 
of [Ket09], by bisecting in r and solving a sequence of linear programming feasibility problems. 
The solution of this optimization problem would be the formula for use in the next step. We do 
not pursue that approach here. Instead, we derive families of formulae that can be applied based 
on the sequence of previous step sizes. 

2.3 Two optimal methods 

For convenience, here we list the methods most likely to be of interest for practical application. 
These (and other methods) are derived and analyzed in the rest of the paper. Recall that p„ has 
been defined in (14), and we assume n > k. The definition of the O quantities (with dependence 
on n suppressed) is given in (4)-(5). 

2.3.1 SSPMSV32 

Our three-step, second-order method is 
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with SSP step size restriction 


hn-2 + h„-l 


/ I 7 I ' 

hn-2 + h„-i + jin 


If the step size is constant, this is equivalent to the known optimal second-order three-step SSP 
method. 


2.3.2 SSPMSV43 

Our four-step, third-order method is 

(03+ 1)2(03-2) 


o| 


W77-I + 


Oq 


3O3 -|- 2 


o| 


Wn-4 + 


Q 3-2 

+ 1 ) 


hn f [U-n—l)^ -\- 


with SSP step size restriction 


h„ < 


3 Q 3 -|- 2 


h„-j 


hn-j] 3“ 


hn f (Wfj— 4 ) 




If the step size is constant, this is equivalent to the known optimal third-order four-step SSP 
method. 


3 Existence and construction of optimal SSP formulae 

In this section we consider the set of formulae satisfying (17) for fixed order p, number of steps 
k and some step-size sequence Clj. It is natural to ask whether any such formula exists, what 
the supremum of achievable r values is (i.e., the optimal SSP coefficient C), and whether that 
supremum is attained by some formula. Here we give answers for certain classes. 

In Section 3.1 we discuss how large an SSP coefficient can be, and prove the existence of a 
formula with the maximum SSP coefficient. In Sections 3.2 and 3.3 we construct some practical 
optimal formulae of order 2 and 3, while the existence of higher-order formulae is established in 
Section 3.4. The theorems of the present section are proved in Section 8. Our theorems are based 
on [NK16] by extending the corresponding results of that paper to the variable step-size case. 
The basic tools in [NK16] include Parkas' lemma, the duality principle, and the strong duality 
theorem of linear programming [Sch98]. 

3.1 Upper bound on the SSP coefficient and existence of an optimal formula 

In the fixed-step-size case, the classical upper bound C < on the SSP coefficient C for a k-step 
explicit linear multistep formula of order p with k > p was proved in [Len89] together with the 
existence of optimal methods. 
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Theorem 2. Suppose that some time-step ratios coj are given. Then the SSP coefficient for a k-step explicit 
linear multistep formula with order of accuracy p >2is bounded by 


C{co,S,f) < 




ifClk < p 
ifClk > p. 


(18) 


Moreover, suppose there exists a k-step explicit linear multistep formula of order p >2 with positive SSP 
coefficient. Then there is a k-step formula of order p whose SSP coefficient is equal to the optimal one. 


3.2 Second-order formulae 


The bound in Theorem 2 is sharp for p = 2, as the following result shows. 


Theorem 3 (Optimal second-order formulae with k >2 steps). Suppose that some time-step ratios 
LOj are given. Then there exists a second order linear multistep formula with k steps and with positive SSP 
coefficient if and only ifClk > 2. In this case, the optimal formula is 


Uyi 


^k-1 ^ 



O.k-1 — 1 


f {Pn—1 


1 

mi 


^n—k/ 


(19) 


and has SSP coefficient 

C{w,d,l5) = ^^. ( 20 ) 

ilk — I 

Remark 3. If the step size is fixed, method (19) with k>3is equivalent to the optimal k-step, second-order 
SSP method given in [GKSll, Section 8.2.1]: 

with SSP coefficient C = {k — 2)/{k — \). 


3.3 Third-order formulae 

Compared to the family of second-order formulae above, the optimal third-order formulae have 
a more complicated structure. Although we will eventually focus on two relatively simple third- 
order formulae (corresponding to k = A and k = 5), we present complete results in order to 
give a flavor of what may happen in the search for optimal formulae. The following Theorem 4 
characterizes optimal third-order linear multistep formulae and their SSP coefficients, again for 
arbitrary step-size ratios cOj. The theorem also provides an efficient way to find these optimal 
A:-step formulae, since the sets of non-zero formula coefficients, denoted hy Af = N{co,5,f>), are 
explicitly described. First we define 2k —2 quantities and sets that will appear in Theorem 4. 


8 










• For j = 0 : 


TV — T' *^0 := 

Llk — i 

• For 1 < j < k — 2: 

/Q^-Oy-S 2 , 1 \ 

Tj := max —---— + —--— , 

— Clj — 1 coj Qjt — Oy—1 J 

and either Sj := or Sj := {Jy,,6y_i,j6y} or Sj := {<Jy,/3y_i,/3y,/3jt_i}, depending 

on whether the first expression is greater, or the second expression is greater, or the two 
expressions are equal in the above max(. 


• For i = k — V. 


rk-i ■■= 


+ 


(^k-l ^k — 


'k-2 


Sk-1 ■= {h-1, fik-2, fik-l}- 


• The rj.+y quantities for 0 < j < k — 3 are defined below, and 5,+y := {/3y,/3,+i,/3,_i}. For 
any 0 < j < k — 3 we set 

Pk+j{^) •= ~ + ^;+i) + 2 (Ay + Ay+j + l) x — 6, (22) 


where 


/\,r, .— fT/, fT,) 


If 


Ay +1 — (Ay + l) Ay_|_i + 3Ay > 0 
then P/c+y (•) has a unique real root. We define 


or 


Ay < 5 


2V6, 


rk+j ■■= 


the real root of Pk+j 
+00 


if (24) holds 

if (24) does not hold. 


(23) 

(24) 


(25) 


Theorem 4 (Optimal third-order formulae with k >2 steps). Let time-step ratios cOj he given. Then 
the inequality Ojt > 3 is necessary and sufficient for the existence of a third-order, k-step explicit linear 
multistep formula with positive SSP coefficient. For > 3, the optimal SSP coefficient is 


Cico,5,B) = min r,-, 

0<j<2k-3 ‘ 


and the set of non-zero coefficients of an optimal SSP formula satisfies J\f C Si, where the index £ G 
{0,1,... ,2k — 3} is determined hv the relation re = min r,-. If the index £ where the minimum is 

^ 0<7<2Jc-3 ‘ ^ 


attained is not unique and we have min r,- = re, 

0 <;< 2 Jc -3 ‘ ^ 


= ... = re, then AC C 5^ n ... n Se. 
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Remark 4. Let us highlight some differences between the fixed-step-size and the VSS cases concerning the 
sets of non-zero formula coefficients. 

1. The pattern of non-zero coefficients for optimal third- or higher-order formulae can be different 
from that of their fixed-step-size counterparts (this phenomenon does not occur in the class of optimal 
second-order formulae). A simple example is provided by the optimal 3rd-order, 5-step formula with 

( 7 11 \ 

(Oo,..., Os) := 

where J\f — {for fi, fi} (the coefficient pattern being similar to the case of the optimal fixed-step-size 
3rd-order, 6-step method). 

2. If the index I with C{cv,S,(3) = r^ is not unique, then the optimal formula has less than 3 non-zero 
coefficients in the general case. 

3. If the index I with C{co,S,(3) = r^ satisfies 1 < £ < k — 2 and the expressions in the max(...) 

are equal, then the optimal formula is generally not unique and has more than 3 non-zero coefficients. For 
example, with ooj := 5 for 1 < j < k — 2, cok-i ■= 4 and cok := 1, we have a one-parameter family of 
optimal methods, and J\f C fk-i} 

4-2 := “ ^fk-i)> fk -3 ■= “ 5) and := ^(20 — fk-i) 

for any 5/16 < fk-i < 2/5. However, for any fixed co it can he shown that there is an optimal formula 
that has at most p non-zero coefficients just as in the fixed-step-size case [Len89]. 

The optimal fixed-step-size SSP method of order 3 and k = A or k = 5 steps has non-zero 
coefficients {4/j6o/jSjc-i} [GKSll, Section 8.2.2], In the rest of this section we consider formulae 
with 4 or 5 steps that generalize the corresponding fixed-step-size methods. A continuity argu¬ 
ment shows that the set of non-zero coefficients is preserved if the step sizes are perturbed by a 
small enough amount. Hence we solve the VSS order conditions (17) with p = 3 for r, Sq, fo and 
fk-T by using Qjt-i > 0/ we obtain that the unique solution is 


r = 


4 = 


^k-l — 2 

Ci.k-1 

^{^k-l + 1 ) “ ^k-1 

03 ' 


fo = 


Clk-i +1 
02 


fk-1 = 


(Ofc_l -|- 1)^ 

ni-i 


The resulting VSS formula reads 

_ -h l)2(0^_i — 2) 

“ 03 


I (Gfc-1 + 1)^ 7 ^ I 

Un-1 H-7^2- f (u„_i)-|- 


' ^n—k ^9 hfij(U„_k)- 


03 


^k-1 


(26a) 

(26b) 


(27) 
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Proposition 1. For 


(28) 


2 < < 2(1 + V2) « 4.828, 

the SSP coefficient of (27) is optimal, and is equal to (26a). 

Proof. By using (26), formula (27) takes the form 


U„ = rfk-lUn-l + fk-lK f{u„-l) + (rfo + ^o)Wn-jt + f{Un-k)- 

By definition (see [GKSll, Chapter 8]), its SSP coefficient is given by 

); ■ f ^l^k-1 _ . f^k-l~'^ 3nj;_i+2 

C [CO, 0, p) = mm ( -, -—^- | = mm 


fk-i ' fo 


from which we see that 


ni-_i —2 


' ^k-l{^k-l +'^) J ' 


for 2 < Clk-i < 2(1 + V2) 


C(co,S,B) = { ,, 

1 n.S‘nt.,1) for n,_, > 2(1 + V^). 


(29) 


But Theorem 2 says that the SSP coefficient of any multistep formula with p = 3 can be at most 
so for 2 < < 2(1 + \/2) the SSP coefficient of (27) is optimal. □ 

The natural requirement (28) also justifies our choice for k: in the fixed-step-size case we have 
— k — 1, and 2 < — 1 < 2(1 -|- -\/2) holds if and only if /c = 4 or /c = 5. 

Remark 5. For k G {4,5}, we have the following strengthening of Proposition 1: the VSS formula (27) 
is optimal if and only if (28) holds. To see this, it is enough to show that (27) is not optimal for 

n,t-i >2(1 + V2). (30) 


Indeed, by fixing any Qjc-i >2(1 + ^), one checks by direct computation that 


30^-1 -|- 2 ^ 

+ 1 ) ^ 


a = 0,1,...,2k-3). 


(31) 


But the SSP coefficient of (27) is given by the left-hand side of (31) according to (29), and the optimal 
SSP coefficient for third-order formulae is min rj according to Theorem 4. Hence the SSP coefficient 

of (27) is not optimal when (30) holds. 


Remark 6. One could develop optimal third-order explicit SSP formulae for k > 5as well. However, their 
structure, as indicated by Theorem 4, would be more complicated, and the analysis performed in Section 
4.2 would become increasingly involved. 
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3.4 Higher-order formulae 

The last theorem in this section reveals that arbitrarily high-order VSS SSP explicit linear multi- 
step formulae exist, though they may require a large number of steps. 

Theorem 5. Let Ki,K 2 > 1 be arbitrary and let ip and k > p^K\K2/2 be arbitrary positive integers. 
Suppose that coj are given and that 

l/Xi < cOj < K 2 for all 1 < j < k. (32) 

Then there exists a k-step formula of order p with C{co,S,f) > 0. 

4 Step-size selection and asymptotic behavior of the step sizes 

To fully specify a method, we need not only a set of multistep formulae but also a prescription 
for the step size. When using a one-step SSP method to integrate a hyperbolic PDE, usually one 
chooses the step size h„ := yC where 7 is a safety factor slightly less than unity. For 

SSP multistep methods, the choice of step size is more complicated. First, multiple previous steps 
must be taken into account when determining an appropriate /ife/ as already noted. But more 
significantly, the SSP coefficient depends on the method coefficients, while the method coef¬ 
ficients depend on the choice of h„. These coupled relations result in a step-size restriction that 
is a nonlinear function of recent step sizes. In this section we propose a greedy step-size selec¬ 
tion algorithm and investigate the d}mamics of the resulting step-size recursion for the formulae 
derived in the previous section. 

Besides the step-size algorithms themselves, our main result will be that the step size remains 
bounded away from zero, so the computation is guaranteed to terminate. Because the step-size 
sequence is given by a recursion involving Izfe/ we will at times require assumptions on hj;^: 

For all n we have p ^ hpEiUn) < for some p^ G (0,oo). (33) 

For all n we have (5 ee < ^ gome prescribed value ()fe S (0,1]. (34) 

“FE(^n+lj (IFE 

Assumption (33) states that the forward Fuler permissible step size remains bounded and is also 
bounded away from zero. For stable hyperbolic PDF discretizations, this is very reasonable since 
it means that the maximum wave speed remains finite and non-zero. Assumption (34) states 
that the forward Fuler step size changes little over a single numerical time step. Typically, this is 
reasonable since it is a necessary condition for the numerical solution to be accurate. It can easily 
be checked a posteriori. 
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4.1 Second-order methods 


Let us first analyze the three-step, second-order method in detail. 

Set k = 3m the second-order formula (19), and suppose that hj>0 has already been defined 
for 1 < i < n — 1 with some n >3. The SSP step-size restriction (12) is implicit, since C„ depends 
on h„. By (20) we have 

^ _ ^l,n 9 “ ^2,n 1 _ 2 9 “ 1 

<^l,n + <^2,n + ^w-1 

Solving for //„ in (12) gives 

j ^ hn—2 hyi—i 
kn ^ ", ' J j ■ }ln- 

hn-2 + kn-l + fin 

It is natural to take the largest allowed step size, i.e., to define 


, hn-2 + hn-l 

kn ■= -j - —j -^- 

«)!-2 + "n-1 + Mn 


■ fin- 


(35) 


For the general k-step, second-order formula (19), the same analysis leads to the following choice 
of step size, which guarantees monotonicity: 


fi . 


hn - — 


fl . 




(36) 


Note that this definition automatically ensures > 1/ arid hence C„ > 0 for any ]^n ^ 0 * 


4.1.1 Asymptotic behavior of the step size 

Since (36) is a nonlinear recurrence, one might wonder if the step size could be driven to zero, 
preventing termination of the integration process. The following theorem shows that, under 
some natural assumptions, this cannot happen. 

Theorem 6. Consider the solution of (1) by the second-order formula (19) with some k>3. Let the initial 
k — 1 step sizes be positive and let the subsequent step sizes h„ be chosen according to (36). Assume that 
(33) holds with some constants Then the step-size sequence h„ satisfies 

k — 2 k — 2 

- - u~ < liminf/z„ < limsup//„ < -- u^. (37) 

k-1 n^+oo k-1 

As a special case, ifh-pEiun) is constant, then 

k — 2 

kn -—-hfE {n -|-oo). 

Remark 7. The asymptotic step size /ipE given above is precisely the allowable step size for the fixed- 
step-size SSP method ofk steps. 
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The proof of Theorem 6 is given in Section 9. 


Remark 8 . Our greedy step-size selection (36) for second-order methods is optimal in the following sense. 
Let us assume that there is another step-size sequence, say h:^^, with the following properties: 

• the corresponding starting values are equal, that is, hj = hj > 0 for j = 1,2,... ,k — 1; 

• the fi„ quantities in (14) corresponding to the sequences h„ and are all equal to a fixed common 
constant y > 0; 

• h^ satisfies (36) with inequality, that is. 


K < 


yk-l 

^;=1 “ n -/ 


Lp; h-_.)+,, 


■y. 


Then—as a straightforward modification of the proof of (78) in Section 9.2 shows—we have < //„ for 
all n > 1. 


4.2 Third-order methods 


In this section we give our step-size selection algorithm for the 3rd-order 4-step and 5-step SSP 
formulae (27) by following the same approach as in Section 4.1. 

By using the optimal SSP coefficient C„ = — given in Proposition 1 we see that h„ = C„y„ 


in ( 12 ) if and only if 


This relation also yields that 


hfi — 


w'c-i h ■ 


YLj—l ^n—j ) T 




■ Pn- 


k-1 


^k—l,n — 2 -|- ^n—jr 

Pn 

SO > 2 is guaranteed by pn > 0. Therefore, (28) is equivalent to 

k-1 

y~! hn—j — Pfi. 

M 


(38) 


The definition of h„ in Theorem 7 below is based on these considerations. The assumptions of 
the theorem on the starting values and on the problem (involving the boundedness of the //fe 
quantities and that of their ratios) are constructed to ensure (38). As a conclusion. Theorem 7 
uses the maximum allowable SSP step size (12) together with the optimal SSP coefficient C„ > 0. 

Theorem 7. Consider the solution of (1) by the third-order formula (27) with k = A or k = 5. Let the 
initial k—1 step sizes be positive and let the subsequent step sizes h„ be chosen according to 


4 ^ 7=1 "-n-j 


hn : — 


yk-ly . 


2pn 


Pn- 


(39) 
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Assume that (33)-(34) and the condition 


0 < hj < Q ■ h^E{uj) for i = 


hold with 


{q,Qfe) = 


(to'to) /or/: = 4, 

/ _ 57 _ 962 '' 1, _ c 

UOO'IOOO/ A. — J. 


Then the step-size sequence hn satisfies 


k-3 


V 


< lim inf h„ < lim sup //„ < 

«—> + 00 


k-3 


F- 


(40) 

(41) 


(42) 


As a special case, ifhEE{un) is constant, then 

k — 3 

h„ -7- — -h^E {n —>• +oo). 

The proof of this theorem is given in Section 9. From the proof we will see that it is possible 
to slightly adjust the simple ((), (5 fe) values given in (41) based on properties of the problem to 
be solved; see Figure 5. 

Remark 9. Our greedy step-size selection (39) for the third-order methods in the above theorem is optimal 
in the same sense as it is described in Remark 8. 


5 Stability and convergence 

The conditions for a method to have positive SSP coefficient are closely related to sufficient 
conditions to ensure stability and convergence. 

Recall that a linear multistep method is said to be zero-stable if it produces a bounded se¬ 
quence of values when applied to the initial value problem 

u'(t) = 0 M(fo) = Uq t S [io/to + T] (43) 

(see, e.g. [HNW93, Section III.5]). The following result shows that variable step-size SSP methods 
are zero-stable as long as the step-size restriction for monotonicity is observed. 

Theorem 8. Let 0 Lj^„ be the coefficients of a variable step-size linear multistep method (3) and suppose that 
^j,n > ^fof n,]. Then the method is zero-stable. 

Proof Application of (3) to (43) yields the recursion 

k-l 

^ ^j,n^n—k+j- 
;=0 

Since the coefficients a.j^„ are non-negative and sum to one, each solution value is a convex 
combination of previous solution values. Therefore the sequence is bounded. □ 
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Corollary 1. Let an SSP variable step-size LMM be given and suppose that the step sizes are chosen so 
that C„ > Q for each n. Then the method is zero-stable. 

In order to prove convergence, we must also bound the local error by bounding the ratio of 
successive step sizes and bounding the method coefficients. 

Lemma 1. For method (19) with k>3or method (27) with k = AorS, let the step sizes be chosen so that 
C„ > 0. Then there exists a constant A such that 

0 < „ < A for all n. 

Proof For any method, C„ > 0 implies that > 0 and 0 < aj)„ < 1. For the second-order 
methods, > 0 implies njc-i,n > 1/ which implies /S^-i ,n < 2. For the third-order methods, 

Cn > 0 implies > 2, which implies fk-i,n <9/4 and /3o,« <3/4. □ 

We recall the following result from [HNW93, Section 111.5, Theorem 5.8]. 

Theorem 9. Let a variable step-size LMM be applied to an initial value problem with a given f and on a 
time interval [to, fo -|- Tj. Let h = max„ h„ be the largest step in the step sequence. Assume that 

1. the method is stable, of order p, and the coefficients (Xj^„, fij^„ are bounded uniformly ash ^ 0; 

2. the starting values satisfy \\u{tj) — Uj\\ = Olh^) where ho is a bound on the starting step sizes; 

3. hn/hn-\ < r] where rj is independent ofn and h. 

Then the method is convergent of order p; i.e., \\u{t„) — u„\\ = 0{hP) for all t„ G [fo,fo + T]. 

The methods described in Section 4 satisfy the conditions of Theorem 9, and are thus conver¬ 
gent, as shown in the following theorem. 

Theorem 10. Under the assumptions of Theorem 6, our methods defined in Section 4.1 and with \\u{tj) — 

UjW = 0{hf) for j < k are convergent of order 2. Similarly, under the assumptions of Theorem 7, our 
methods defined in Section 4.2 and with \\u{tj) — Uj\\ — 0{h^) for j <k are convergent of order 3. 

Proof. It is sufficient to show that assumptions 1 and 3 of Theorem 9 are fulfilled. Our construc¬ 
tion in Sections 4.1 and 4.2 guarantees > 0. Therefore Lemma 1 applies, so the coefficients 
f>j,n are uniformly bounded and non-negative. Thus Theorem 8 applies, so the methods are 
stable. Furthermore, the condition h„/h„_i < rj < -|-oo is implied by (37) or (42), hence Theorem 
9 is applicable. □ 

6 Numerical examples 

In this section we investigate the performance of the proposed methods by performing numer¬ 
ical tests. In Section 6.4, the accuracy of our methods is verified by a convergence test on a 
linear equation with time-varying advection velocity. In Sections 6.5-6.8 we apply the methods 
SSPMSV32 and SSPMSV43 to nonlinear hyperbolic conservation laws in one and two dimensions. 

All code used to generate the numerical results in this work is available at https: //github. com/numerical-m 
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6.1 Efficiency 

Let N denote the number of steps (excluding the starting steps) used to march with a fc-step 
variable step-size method from time to to a final time to + T. If the same method were used 
with a fixed step size, the number of steps required would be at least N' = {T — jj-Zl 
where limin is the smallest step size used by the variable step-size method. Thus the reduction in 
computational cost by allowing a variable step size is given by 


■ N' /Zavg' 

where havg := (T — hj)/N is the mean step size used by the variable step-size method. 

6.2 Spatial discretization 

In space, we use the wave-propagation semi-discretizations described in [KPL13]. For the 
second-order temporal schemes, we use a spatial reconstruction based on the total-variation- 
diminishing (TVD) monotonized-central-difference (MC) limiter [Lee77]; for the third-order tem¬ 
poral schemes, we use fifth-order WENO reconstruction. The only exception is the Woodward- 
Colella problem, where we use the TVD spatial discretization for both second- and third-order 
methods. The second-order spatial semi-discretization is provably TVD (for scalar, ID problems) 
under forward Euler integration with a CEL number of vpE = 1/2. We use this value in the 
step-size selection algorithm for all methods and all problems. 

6.3 Time-stepping algorithm 

The complete time-stepping algorithm is given in Algorithms 1 and 2. We denote by the 
CEL number at step n. By default we take the initial step size hi = 0.1. In order to ensure 
monotonicity of the starting steps, we use the two-stage, second-order SSP Runge-Kutta method 
to compute them, with step size 

h„ := yCohniun-i), n = l,...,k-l, 

where Co = 1 is the SSP coefficient of the Runge-Kutta method and 7 = 0.9 is a safety factor (cf. 
Section 4). 

For our third-order methods we can check conditions (34) and (40) only a posteriori. If con¬ 
dition (34) or (40) is violated, then the computed solution is discarded and the current step is 
repeated with a smaller step size as described in Algorithm 2. These step-size reductions can be 
repeated if necessary, but in our computations the first reduction was always sufficient so that 
the new step was accepted. For the numerical examples of this section, the step size and CEL 
number corresponding to the starting methods are shown on the right of Figures 1 and 2. 
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Algorithm 1 Second-order time-stepping algorithm 

Initialization: Choose an initial step size hi and compute the forward Euler step size //fe(mo). 

Starting procedure: 

for n = 1,2,..., A: — 1 do 

Compute u„ using the two-stage, second-order SSP Runge-Kutta method with 
step size h„. 

Find the CEL number v„ and compute 
Set h ^ -fCohi^{u„). 
if v„ > Co i/pE then 

Set h„ h and repeat the step, 
else if n < A: — 1 then 
Set hyij^i i — h. 
else 

Set Az/c based on (36). 

end if 
end for 

Main method: 
for n = k,. .. ,N do 

Compute Un using the method (19) with step size /z„. 

Compute from (36). 

end for 


Remark 10. There is no need to check the CFL condition i/„ < Co Vfe for n > k (i.e. when the LMM 
method is used) since it is automatically satisfied by the step size selection (12). 

Condition (34) is violated only when the maximum wave speed changes dramatically between 
consecutive steps, which may suggest insufficient temporal resolution of the problem. We have 
found that, for the problems considered herein, omitting the enforcement of condition (40) never 
seems to change the computed solution in a significant way. Since these two conditions were 
introduced only as technical assumptions for some of our theoretical results, they could perhaps 
be omitted in a practical implementation. 

6.4 Convergence test 

We consider the linear advection problem 

ut+ + ^sm(27rf)^ 
u{x,0) = sm{2nx), x 


Wx = 0, (44) 

G [0,1], 
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Algorithm 2 Third-order time-stepping algorithm 

Initialization: Choose an initial step size hi and compute the forward Euler step size //fe(mo)- 

Starting procedure: 

for n = 1,2,..., A: — 1 do 

Compute u„ using the two-stage, second-order SSP Runge-Kutta method with 
step size h„. 

Find the CEL number v„ and compute //fe(wm). 
if condition (34) or (40) does not hold then set 


hfi 


hni2 

7Co (?Efe(Wh) 


if at least (34) is violated; 
if (40) is violated. 


and repeat the step, 
else 

SetE ^ 7 C 0 huiun). 

end if 

if Vn > Co vfe then 

Set h ,2 <— h and repeat the step, 
else if n < A: — 1 then 
Set //„+i ^ h. 

else 

Set hk based on (39). 

end if 
end for 

Main method: 
for n = k,... ,N do 

Compute u„ using the method (27) with step size h„. 
if condition (34) is violated then 
Set h„ hn/2 and repeat step. 

else 

Compute hn+i from (39). 

end if 
end for 


with periodic boundary conditions. We use a spatial step size Ax = 2^^, d = 6 ,..., 11, and 
compute the solution at a final time t = 5 (i.e., after 10 cycles). Table 1 shows the Li-norm of the 
error at the final time. The solution is computed by using second-order methods with A: = 3,4 
steps and third-order methods with A: = 4,5 steps. All methods attain the expected order. 
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N 

SSPMSV32 


SSPMSV42 


SSPMSV43 


SSPMSV53 


128 

1.50 X 10-2 


1.83 X 10-2 


9.20 X 10-3 


6.08 X 10-3 


256 

4.30 X 10-3 

1.80 

5.34 X 10-3 

1.78 

1.30 X 10-3 

2.82 

8.10 X 10-3 

2.91 

512 

1.15 X 10-3 

1.90 

1.44 X 10-3 

1.89 

1.68 X 10-2 

2.95 

1.04 X 10-3 

2.96 

1024 

3.01 X 10-4 

1.93 

3.81 X 10-4 

1.92 

2.13 X 10-3 

2.98 

1.32 X 10-2 

2.98 

2048 

7.74 X 10-3 

1.96 

9.84 X 10-3 

1.95 

2.67X 10-‘^ 

2.99 

1.66 X 10-3 

2.99 


Table 1: Li-norms of the error at final time t = 5 for the variable-coefficient advection prob¬ 
lem (44) for second- and third-order LMMs. For each method, the second column denotes the 
convergence order. The number of spatial points is indicated by N. 


6.5 Burgers' equation 

We consider the inviscid Burgers' initial value problem 


Ut + UUx = 0, 

1 

u{x,0) = - -I-sin(27rx), x G [0,1], 

with periodic boundary conditions and 256 grid cells. Figure la shows the evolution of the 
step size, up to a final time t = 0.8. The step size is constant until the shock forms, and then 
increases. The dashed curves in Figure la indicate how the CFL number varies with time for 
each method. Since t/„ = C„Vfe and vpE = 1/2, the CFL number after the first few steps is 1/4 
and 1/6, for the second- and third-order LMMs, respectively. Note that these are the theoretical 
maximum values for which the solution remains TVD. For the third-order method coupled with 
WENO discretization the TVD-norm of the solution does not increase more than 10^^ over a 
single time step. The sudden decrease of the step size after the first few steps of the simulation 
is due to the switch from the starting (Runge-Kutta) method to the linear multistep method. For 
both methods we have s ~ 0.88. 


6.6 Woodward-Colella blast-wave problem 

Next we consider the one-dimensional Euler equations for inviscid, compressible flow: 

Pt + {pu)x = 0 
{pu)t + {pu^ + p)x = 0 

Et + {u{E -p p))x = 0. 

Here p, u, E and p denote the density, velocity, total energy and pressure, respectively. The fluid 
is an ideal gas, thus the pressure is given by p = p{'Y — l)e, where e is internal energy and 
7 = 1.4. The domain is the unit interval with Ax = 1/512, and we solve the Woodward-Colella 
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blast-wave problem [WC84] with initial conditions 


p{x,0) = l, v{x,0)=0, 


and 


p{x,0) = < 


1000 if 0 < X < 0.1 
0.01 if 0.1 < z < 0.9 . 

100 if 0.9 < X < 1 


The initial conditions consist of two discontinuities. Reflecting boundary conditions are used and 
initial conditions lead to strong shock waves, and rarefactions that eventually interact with each 
other. The second-order TVD semi-discretization is used with both the second- and third-order 
methods and the solution is computed at time t = 0.04, after the two shock waves have collided. 

Figure lb shows how the step size evolves over time. The maximum wave speed decreases 
as the shock waves approach each other, then increases when they collide, and finally decreases 
when the shocks move apart. The step size exhibits exactly the opposite behavior since it is 
inversely proportional to the maximum wave speed. As before, the CFL number remains 
close to 1 /4 and 1/6 for the second- and third-order LMMs, respectively. Note that these are the 
theoretical maximum values for which the characteristic variables remain TVD. For this problem 
s 0.76. 


6.7 Two-dimensional shallow-water flow 

Consider the two-dimensional shallow-water equations 

ht + {hu)x + {hv)y = 0 
{hu)t + {hu^ + ]^gh^)x + {huv)y = 0 

{hv)t + ihuv)x + {hv^ + = 0, 

where h is the depth, {u,v) the velocity vector and hu, hv the momenta in each direction. The 
gravitational constant g is equal to unity. The domain is a square [—2.5,2.5] x [—2.5,2.5] with 
250 grid cells in each direction. The initial condition consists of an uplifted, inward-flowing 
cylindrically symmetric perturbation given by 

/i(x,y,0) = l+g{x,y), u{x,y,0) = -xg{x,y) and v{x,y,0) = -yg{x,y), 

where g{x,y) := e . We apply reflecting boundary conditions at the top and right 

boundary, whereas the bottom and left boundaries are outflow. The wave is initially directed to¬ 
wards the center of the domain, and a large peak forms as it converges there, before subsequently 
expanding outward. We compute the solution up to time t = 2.5, after the solution has been re¬ 
flected from the top and right boundaries. Figure 2a shows how the step size varies in time. It 
decreases as the initial profile propagates towards the center of the domain since the maximum 
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(b) Woodward-Colella blast-wave problem 


Figure 1: Evolution of step size for one-dimensional Burgers' and Euler equations with second- 
and third-order linear multistep methods (solid lines). The dashed curves indicate the CFL 
number for each method. Figures on the right show a close-up view of the step size and CFL 
number at early times. The black marks indicate the starting Runge-Kutta method's step sizes. 


wave speed increases. As the solution expands outwards, higher step sizes are allowed. Note 
that at time t k, 2 the wave hits the boundaries, so a small decrease in the allowed step size is 
observed. Again, the dashed curves indicate the variation of the CFL number. Now we have the 
ratio s ~ 0.62 for both second- and third-order methods. 

6.8 Shock-bubble interaction 

Finally, we consider the Euler equations of compressible fluid dynamics in cylindrical coordi¬ 
nates. Originally, the problem is three-dimensional but can been reduced to two dimensions by 
using cylindrical symmetry. The resulting system of equations is 
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pt + {pu)z + {pv)r - -y 
{pu)t + {pU^ + p)z + (|OUZ?)r = 

(pz;)t + + {pv^ + p), = -^ 

{pE)t + {u{pE + p))z + iv{pE + p)), = 

Here, p is the density, p is the pressure, E is the total energy, while u and v are the z- and 
r-components of the velocity The z- and r-axis are parallel and perpendicular, respectively, to 
the axis of symmetry The problem involves a planar shock wave traveling in the z-direction 
that impacts a spherical low-density bubble. The initial conditions are taken from [Ket-i-12]. 
We consider a cylindrical domain [0,2] x [0,0.5] using a 640 x 160 grid and impose reflecting 
boundary conditions at the bottom of the domain and outflow conditions at the top and right 
boundaries. 

Since the fluid is a polytropic ideal gas, the internal energy e is proportional to the temper¬ 
ature. Initially the bubble is at rest, hence there is no difference in pressure between the bubble 
and the surrounding air. Therefore, the temperature inside the bubble is high resulting in large 
sound speed c = \/j(7 — l)e, where 7 = 1.4. Consequently, this results in reduced step sizes 
right after the shock wave hits the bubble at f ~ 0.05, as shown in Figure 2b. The efficiency ratio 
s is about 0.82. 

7 Conclusions 

The methods presented in this paper are the first adaptive multistep methods that provably 
preserve any convex monotonicity property satisfied under explicit Euler integration. We have 
provided a step-size selection algorithm—yielding an optimal step-size sequence hn —that strictly 
enforces this property while ensuring that the step size is bounded away from zero. The methods 
are proved to converge at the expected rate for any ODE system satisfying the forward Euler 
condition. 

As suggested by Theorems 6 and 7, for all tests the CEL number remains close to IrfwE/ 
where k and p are the steps and order of the method, respectively. The numerical results verify 
that the step size is approximately inversely proportional to the maximum wave speed and the 
proposed step-size strategy successfully chooses step sizes closely related to the maximum al¬ 
lowed CEL number. In practice, we expect that conditions (34) and (40), which were introduced 
only as technical assumptions for some of our theoretical results, could usually be omitted with¬ 
out impacting the solution. 

The methods presented herein are of orders two and three. We have proved the existence 
of methods of arbitrary order. However, the optimal methods for orders at least 3 seem to 
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Figure 2: Evolution of step size for two-dimensional shallow-water and Euler equations with 
second- and third-order linear multistep methods (solid lines). The dashed curves indicate the 
CFL number for each method. Figures on the right show a close-up view of the step size and 
CFL number at early times. The black marks indicate the starting Runge-Kutta method's step 
sizes. 


have a complicated structure. It would be useful to develop adaptive (and possibly suboptimal) 
multistep SSP methods of order higher than three having a simple structure. 


8 The proofs of the theorems in Section 3 

We present three lemmas. Lemma 2 below. Lemma 3 in Section 8.1 and Lemma 4 in Section 8.2, 
that are used in the proofs of Theorems 2-5 in Sections 8.1-8.4. Lemmas 2-4 are straightforward 
generalizations of the corresponding results of [NK16] for variable step-size formulae, hence their 


24 



















































proofs are omitted here. 

Lemma 2. Let r > 0 be arbitrary and let p and n > kbe arbitrary positive integers. Suppose that some 
time-step ratios C 0 j^„ > 0 are given. Then the following two statements are equivalent. 

(i) For all formulae with k steps and order of accuracy p we have Cn{oo, 5, f>) < r. 

(ii) There exists a non-zero real univariate polynomial q„ of degree at most p, that satisfies the conditions 


q„{n,j^„) >0 for all 0 < j <k — 1, (45a) 

q'n{Clj^n) + >0 for all 0 < j < k — 1, (45b) 

hn{k^k,n)=0. (45c) 

Furthermore, if statement (z) holds, then the polynomial q„ can be chosen to satisfy the following condi¬ 
tions. 

1. The degree of q„ is exactly p. 

2. The real parts of all roots ofq„ except for interval [Oo^n, 

3. The set of all real roots ofq„ is a subset of {Oo^n, ..., Oj. 

4. For even p, is a root of q„ with odd multiplicity. For odd p, if is a root of q„, then its 
multiplicity is even. 

5. The multiplicity of the root is one. For 1 < j < k — 1, if Cij^n is a root of q„, then its 
multiplicity is 2. 

6. The polynomial q„ is non-negative on the interval [Oo^n, 

The following observation will be useful in the proofs. If q„{Cij^„) = 0, then the inequality 
(45b) for this index j simplifies to q'„{aj^„) > 0. Otherwise, if q„{Clj^„) 0, then the inequality 

(45b) for this index; can be written by using the logarithmic derivative of qn as 

^ 1 ^ 

qn{F\jn) 4.^,; 

where the „ numbers are the complex (including real) roots of qn, and qn has been chosen such 
that its degree is exactly p. 
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8.1 The proof of Theorem 2 

The following lemma will be used in the second step of the proof of Theorem 2. 

Lemma 3. Let r > 0 be arbitrary and let p and n > kbe arbitrary positive integers. Suppose that some 
time-step ratios C 0 j^„ > 0 are given. Then exactly one of the following statements is true. 

(i) There is a formula with k steps, order of accuracy p and C„{co, S, f>) > r. 

(ii) There is a real univariate polynomial q„ of degree at most p, that satisfies the conditions (45a), (45b) 
and 

hn{L^k,n) ^ 0- (46) 


The proof of Theorem 2. Step 1. To prove the bound (18), we first consider the k >2 case, and let 
qn denote the polynomial 

qn{x)-.= x). (47) 

• If O/f > p, this qn satisfies conditions (45) with r := ~ p) / {^k,n ~ 1)/ since for each 

1 < j < k — 1 we have 


qn {Clj n) 


P - 1 _ 1 ^ P - 1 

^j,n ^k,n ^j,n ^k—l,n 


Q 


k,n 


1 

— ^k-l,n 


^k,n P 

^k,n — 1 ' 


and q'n{O. 0 ,n) = q'n{0) > 0. 

• If n.k,n < P/ the conditions (45) hold with r := 0, because for x G = [0, ~ 1] 

we have 

q',fx) > xP-^ ((p - l)nk,n - V {^k,n - 1)) = xP-^ (p - n^^n) > 0. 

For k = \, the same qn as in (47) satisfies (45) with r := 0 and Q;t,« = 1- Then, in each case, (18) is 
an immediate consequence of Lemma 2. 

Step 2. Suppose now 


3 a fc-step explicit linear multistep formula of order p > 2 with Cn{co,5,f>) > 0. (48) 


We prove that there is a fc-step formula of order p whose SSP coefficient is equal to the optimal 
one. 

Let us set Tf := {r > 0 : the statement {ii) of Lemma 3 holds}. Step 1 implies that the largest 
possible Cn{(V,5,f>) for all fc-step formulae of order p and having the given time-step ratios coj is 
finite. So 7 ^ 0 by Lemma 3. We also see that "H is an infinite interval, because if a polynomial 
qn satisfies the conditions (45a), (45b) with some r := p > 0, and (46), then it satisfies the same 
conditions with any r G (p, -|-oo). Thus, with a suitable r* > 0, we have (r*, -|-oo) C "H C [r*, -|-oo). 
Clearly, r* > 0 due to assumption (48) and Lemma 3. We claim that r* ^ TL. 
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Suppose to the contrary that r* El-i. Then there is a polynomial cfn satisfying the conditions 
(45a), (45b) with r := r*, and (46). Now we define := (j,, + \qn{^k,n) \ /2. One easily checks 
that this qn polynomial satisfies the same set of conditions with 

1 + max qn(Clj„) 

1 4“ (O;^,;) I/2 + max (|,;(0;^) 

0<;<Jc-l ' ^ ' 

SO we would get (0, r*) n "H 7 ^ 0, a contradiction. 

Hence T-L = {r* , + 00 ). Therefore, in view of Lemma 3, r* is the optimal SSP coefficient and 
there exists an optimal method with C„{co,5,^) = r*. □ 

8.2 The proof of Theorem 3 

In the following we apply the usual terminology and say that an inequality constraint is binding 
if the inequality holds with equality. The lemma below is used in the proofs of Theorems 3 and 
4. 

Lemma 4. Let n > k be arbitrary positive integers and p > 2. Suppose that the time-step ratios 
coj^n > 0 are given, and that there exists an explicit linear multistep formula with k steps, order of 
accuracy p, and Cn > 0. Let and fj^n denote the coefficients of a formula with the largest SSP 
coefficient Cn{co„, S„, fn) G (0,+ 00 ). Let q„ be a polynomial that satisfies all the conditions of Lemma 2 
with r : = C„{cVn,S„, fn). Then the following statements hold. 

(i) If dj,n f Ofor some 0 < j < k — 1, then q„{Ci.j^n) = 0. 

(ii) If 7^ Ofor some 0 < j < k — 1, then q'„{D,j^n) + rqn{L^i,n) = 0. 

(Hi) This q„ can be chosen so that the total number of binding inequalities in (45a)-(45b) is at least p. 

The proof of Theorem 3. The necessity of the condition > 2 for the existence of a second-order 
formula with positive SSP coefficient is an immediate consequence of Theorem 2. On the other 
hand, we easily see that, up to a positive multiplicative constant, q„{x) = x{Clk^„ — x) is the 
unique pol}momial satisfying (45a), (45c), and Properties 1 and 4 of Lemma 2. If > 2, then 
qn does not satisfy (45b) with r := 0 for j = k — 1, therefore there exists a formula with C„ > 0 
due to Lemma 2. Thus the condition > 2 is also sufficient. 

From now on we assume > 2, and that a formula with the optimal SSP coefficient 

C„{co„,6n, fn) is considered. From the above form of qn we see that qn{Dj,n) = 0 (0 < / < k — 1) 

holds only for j = 0, so the statement (in) of Lemma 4 guarantees that at least one of the 
inequalities in (45b) is binding with r Cn{cOn,S„, fn)- Since q'„(Clo,n) + rqn(Clo,n) = Ptk,n > 0 
and for all 1 < / < k — 2 we have 

_1_ 1 ^1_ 1 _ — 1 
qnif-)j,n) ^j,n ^k,n ^j,n ^k—l,n ^k,n ^k—l,n ^k—l,n 



j e (o,H), 
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the unique binding inequality in (45b) must be the one corresponding to the index j = k — 1. But 

0 implies T ^) / {f^k,n !)• 

Finally, the statements (z) and {ii) of Lemma 4 guarantee that only and can differ 

from zero. Solving the conditions for order two (see (17)) then yields the unique formula stated 
in (19). □ 


8.3 The proof of Theorem 4 


First we present a lemma that will be used in the proof of Theorem 4. Moreover, Statement 2 
of the lemma proves our claim about the unique real root of the polynomial in (22) under 
assumption (24), hence guaranteeing that the quantity in (25) has a proper definition. The 
Ay quantities defined in (23) clearly satisfy the condition 1 < Ay+i < Ay below. 


Lemma 5. Let some 1 < Ay+i < Ay numbers be given, and for any p ^Rlet 


Mp) — 


( P-2 


pAj-2Aj 

yP-^y+i “ 2Ay+i 


p-1 p-3 \ 

pAy-1 pAj-3Aj 

pAy_|_i — 1 pA^^^ — 3A^^^ J 


(49) 


Then the following statements hold. 

1. Every real root of the polynomial P/^+y (appearing in (22)) is positive. 

2. Under condition (24), Pk+j has a unique real root. 

3. If 

3{r,a,b) G satisfying Aj{r) ■ {a,b,\)^ = 0 and a^ 


4fo < 0, 


then 

(24) and Pk+jf) = 0 


hold. 

4. If (51) holds with some r G R, then (50) is true, and the triplet {r,a, b) G R^ is unique. 


(50) 

(51) 


Proof. Throughout the proof we always assume 1 < Ay_|_i < Ay and r G R. First we define some 
auxiliary polynomials—their dependence on Ay+i is suppressed for brevity. 

Pi{r) ■■= (r^ -r)Ay+i -r + 2, 

P 2 (r) := (r2-2r)Ay+i-2r + 6, 

P 3 (r) := (r^ - r^)Aj^.^ +2{r^- 2r)Ay+i - 2r + 6, 

Pi{r) := r^ {f- — 4r + 2) Ay^^ — 4r (r^ — 5r + 3) Ay+i + 2 (r^ — 6r + 6), 
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Ay+1 10 


5 


0 

0 5 10 15 20 

Figure 3: The set defined by (24). 

ft(r) := (- (5 +2^6) r3 + 2 (3 + ^) - 2r) Ay+i+ ^5 + 2^) - 4 (3 + v^) r + 6, 

P(,{r) ■= (r - l)Ay+i - r + 3, 
p7{r) ■= - 4rAy+i + 6, 

Ps{x,a,r) := {r^ — r)x^ + {r — \){ar + r — 3)x + a{2 — r) — r + 3, 
p^[r) := (r - - {r^ - 4r + 3) Ay+i + r-3, 

Pio{r) ■■= {r^ -r)Aj+-i-r + 3, 

Pii(r) := rAj^^ - (r + 3)Ay+i + 1. 

Step 1. Statement 1 follows from the fact that the coefficient of x^ in Pi^^j{x) is positive for 
odd m and negative for even m, so cannot have a non-positive root. 

Step 2. We now prove Statement 2. Since Pk+j is cubic, it has a real root. By taking into 
account Statement 1, we will show that 

3r > 0 : (24), P,+y(r) = 0 and P;+y(r) < 0 (52) 

is false, implying that P^^^ > 0 at every real root of Pk+j, which is possible only if Pj.+y has a 
unique real root. 
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If Pi(r) =0, we obtain 


(53) 


{r,Aj+i) = {2-V2,3 + 2V2) 

by using Pjt+y(r) = 0 also. But then one easily shows that (52) is impossible. 

So we can suppose in the rest of Step 2 that Pi(?') 7 ^ 0. Then from Pjt+y (r) = 0 we get 

A^ = P 2 (r)/(rPi(r)). (54) 

This form of Ay is substituted into (52) and we obtain 

r > 0 and P 3 (r)/Pi(r) > 0 and P 4 (r)/Pi(r) < 0 and (55) 


(P 5 (r)/Pi(r) <0 or P6(r) P 7 (r)/Pi(r) > o) . (56) 

Because all Py polynomials (1 < i < 7) are at most quadratic in Ay+i, one can systematically 
reduce (55)-(56) to a system of univariate polynomial inequalities in r, and verify that (55)-(56) 
has no solution. This finishes the proof of Step 2. 

Step 3. To prove Statement 3, notice first that 

Aj{r)-{a,b,l)^ =0 (57) 

implies 

0 = det Ay(r) = (Ay - l) (Ay+I - 1) (Ay - Ay+i) P,+y(r), (58) 

so we have P)c+y(r) = 0. To show (24), depicted in Figure 3, we separate two cases. 

If r = 1, then we obtain a = —2 from the first component of (57), and b = 4Ay — A? with 

Ay.^_i = 4 — Ay and 2 < Ay < 3 (59) 


from the other two components and from 1 < Ay+i < Ay. So (24) holds. 

We can thus suppose in the rest of Step 3 that r ^ 1. Then from (57) we get 

f, ^ (2-r)fl-r + 3 
r — 1 


(60) 


and 

Ps{Aj,a,r) = 0 = P8(Ay+i,fl,r). (61) 

We again separate two cases. 

If Pi(r) =0, then by using (61) and 1 < Ay+i < Ay we obtain 

r-2-V2, Ay > Ay+i = 3 + 2^2 and fl - - (^Ay + 1 + 1/V 2 ) . (62) 

But it is easy to check that (62), the inequality — 4fo < 0 from (50), and the negation of (24) 

cannot hold simultaneously. So (24) is proved in this case. 
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Therefore, in the rest of Step 3, we can suppose that Pi{r) 7 ^ 0. Then a is expressed from 
Ps{Aj+i, a,r) =0 and we get 

fl = P 9 (?')/Pi(r). (63) 

Moreover, we express Aj from P;t+/(^) = 0 (shown to hold in (58)) as in (54). To finish Step 3, we 
will verify that 

—Ab < 0 from (50), and the negation of (24) (64) 

cannot hold. The substitution of (63) and (54) into (64) yields 

r>0, P6{r)Pio{r)Pn{r) < 0, P 5 {r)/Pi{r) > 0, and P 6 {r)P 7 {r)/Pi{r) < 0. (65) 

Again, these Pj polynomials are at most quadratic in Aj^i, so one can check that (65) has no 
solution. 

Step 4. To prove Statement 4, we first suppose that 

r 7 ^ 1 and Pi(r) 7 ^ 0. ( 66 ) 

Then we have seen in Step 3 that if {a, h) G solves (57), then a and h necessarily have the form 
(63) and (60), respectively; and the uniqueness of r is guaranteed by Step 2. So, under condition 
( 66 ), any triplet solving (57) is unique. On the other hand, direct substitution into (57) shows that 
the above (r, a, b) is indeed a solution—the non-trivial component to check in (57) is the second 
one, which is just [Aj — l) (Ay — Ay_|_i) Pjt+y(r)/Pi(r) = 0 due to Pk+j{r) = 0 from assumption 
(51). Finally, to finish Step 4 in the case when ( 66 ) holds, we notice that for this unique triplet 
{r,a,b) we have — Ab = P6{r)PiQ{r)Pii{r)/ (Pi(r))^. But, with Ay expressed as in (54), one can 
again show that 

P6{^)Pio{^)Pii{^) >0/ 1 < A;+i < arid (24) from assumption (51) 

cannot hold. This means that — Ab < 0 from (50) also holds. 

To prove Statement 4 when r = 1, we recall from Step 3 that in this case the only possible 
triplet is {r,a,h) = (1, —2,4Ay — A?) with (59), which indeed satisfy (50). 

Lastly, we consider the case when Pi (r) = 0. We have seen in Step 2 that this time (53) holds. 
Since Ay+i < Ay, we also have 3 + 2\/2 < Ay. These yield that the unique solution to (57) is 

(r,fl,&) - ( 2 --(Ay + 1 + 1 /V 2 ), (2 + ^)Ay). 

Now (24) implies Ay < 5 + 7/\/2, so — Ab = Aj — 3(2 + \/2) Ay + -\/2 + 3/2 < 0, and the proof 
is complete. □ 

The proof of Theorem 4. Step 1. The necessity of the condition Clj^ „ > 3 for the existence of a third- 
order formula with positive SSP coefficient is again an immediate consequence of Theorem 2. 
To see that this condition is also sufficient, suppose to the contrary that there is no third-order 
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formula with k steps, positive SSP coefficient and > 3, and let us apply Lemma 2 with 
r := 0. First we show that, up to a positive multiplicative constant, qn(x) = — x) is the 

unique polynomial satisfying (45), and Properties 1-2 and 4-5 of Lemma 2. Indeed, (45a)-(45b) 
with i = 0 and r = 0, together with Properties 2 and 5 imply that Qo,n = 0 is a root of qn 
(otherwise the logarithmic derivative of qn at 0 would be negative), and then the multiplicity of 0 
is precisely two due to Properties 1 and 4, so the desired form of qn follows. Now, applying (45b) 
with i = k — 1 and r = 0 to this q„ gives (3 — Oj; „) (0^,^ — 1) >0, contradicting to Cl^^n > 3. This 
contradiction proves the sufficiency of the condition Clj^ „ > 3. 

Step 2. From now on we assume > 3, and that a formula with the optimal SSP coefficient 
r := Cn{c0n,5n,^n) S (0,+co) is considered. Let qn be a polynomial satisfying conditions (45), 
Properties 1-6 of Lemma 2, and the condition {in) of Lemma 4. 

Case 1. In the case when Oj^ n is not the only real root of qn, then, due to Properties 1, 3 
and 5 of Lemma 2, qn is unique (up to a positive multiplicative constant) and takes the form 
qn{x) = {x — Oyg „)^(Q)t,,; — x) with an appropriate index 0 < jo < k — 1. We notice that the 
only binding inequality in (45a) is the one corresponding to j = jo, and the j — jo inequality in 
(45b) is also binding (independently of the value of r). On the other hand, now all roots of qn 
are real, so the logarithmic derivative of qn is strictly decreasing on the intervals [Clo,n, ^jo,n) arid 
{Clj^^n,^k-i,n]- Thus, by {in) of Lemma 4, (45b) contains either two or three binding inequalities, 
and they correspond to 

A. j G {0,k-l} for jo = 0, 

B. j e {jo - l,;o} or j e [jo, fc - 1} for l < /o < fc - 2, 

C. 7 e {k — 2, k — 1} for jo = k — 1. 

In Case A, (45b) is solved with j = k — 1, and we obtain r = {Cij^^n “3)/ {Clk,n “ !)• Due to (/) 
and {ii) of Lemma 4, all formula coefficients except for possibly {5o,n, fio,n, fik-i,n} vanish. 

In Case B, (45b) is solved with j = jo — f and j = k — 1 to get r = 2/(Oyg^n ~ ^jo-i,n) + 
1 / (Ofc^n — ^jo-i,n) and r = {Clk,n — ^io,n “3)/ (nj:,„ — f^jg,n ~ 1 )/ respectively, then the maximum 
is chosen. If these two expressions for r are not equal, we get from {i) and {ii) of Lemma 4 that 
the number of non-zero formula coefficients is at most three; whereas if the two expressions for 
r are equal, the number of non-zero formula coefficients is at most four. 

In Case C, (45b) is solved with j = k — 2 to yield r = 2/(Ojt-i^n — 0 ,^-2,n) + 1/ {^k,n — ^k- 2 ,n)- 
All formula coefficients except for possibly {Sk-i,n, (ik- 2 ,n, l^k-Ut} vanish. 

Therefore, the optimal SSP coefficient is equal to one of the rj quantities (0 < ;■ < k - 1) 
defined in Section 3.3, and the non-zero formula coefficients also have the form stated there. 

Case II. In the case when is the only real root of qn, then, by taking into account (45a), 
and Properties 1 and 5 of Lemma 2, qn (up to a positive multiplicative constant) has the form 

qn{x) {f^k,n x) {{k^k,n x) -|- x) b) 
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with some coefficients a, fo G K satisfying the discriminant condition —Ab < 0. Let us introduce 
the abbreviation Qa := + rq„. 

This time there are no binding inequalities in (45a), so, due to (Hi) of Lemma 4, precisely three 
inequalities are binding in (45b) —there can be no more, since the polynomial Qa is cubic. Let 
0 < /o < ji < — \ denote the three indices corresponding to these binding inequalities, then 

= Q 3 (%,„) = = 0. We observe that the leading coefficient of Qa is -r < 0, 

and Qa(Llj.„) — —b < 0 because of — Ab < 0. As is strictly increasing in j, we get that 
Qa is positive on (-oo,Qyj, „) U and negative on U Hence 

neither ji > jo+ 2 nor j 2 < k — 2 can occur (otherwise < 0 or Qa(n|f_i „) < 0 would 

contradict to (45b)). This shows that the binding inequalities in (45b) are the ones corresponding 
to the index set j G { 70 , jo + l,k — 1} with an appropriate index 0 < jo < k — 3. The previous 
sentence can be rewritten as Aj^ (r) • {a, b, 1)^ = 0, with the matrix Aj^ defined in (49). Since now 
— 4& < 0 as well, so (50) is also satisfied. Hence (24) and Pk+jg {^) = 0 hold by Lemma 5. Due 
to uniqueness we see that r = Tk+jg, with ri+jg defined in (25). 

Thus, the optimal SSP coefficient is equal to one of the quantities (0 < jo < k — 3), and 
we get from (h) of Lemma 4 that all formula coefficients except for possibly {^jg,n, ^jg+i,n, f^k-i,n} 
vanish. 

Step 3. Finally, in order to conclude that the optimal SSP coefficient C„{co„,S„,I3„) is in fact 
the minimum of the rj expressions defined in Section 3.3, we show that C„{co„,3n, fin) < for 
any 0 < j < 2k — 3. 

Indeed, fix any 0 < 70 < A: — 1 and define qn{x) := {x — „)^(njt„ — x). The inequalities 

(45a) and (45c) are trivially satisfied. We verify (45b) with r := ry^ by setting Q 3 ■= q'„ + rq„ and 
distinguishing three cases. 

A. If jo = 0, then one checks that the three roots of Qa are found at {xi,0, with 

Xi = —2Clk^nl {^k,n — 3) < 0. 


B. If 1 < 70 < k — 2 and 

k^k,n k^jo,n 3 2 1 

-^- — - <-1-, 

k^k,n 1 ^jg,n P^k,n 


(67) 


then the three roots of Qa are found at {ny(,_i^„,nyp ,„xa} with Clk-i^„ < X 3 < Clk^n- If 
we have ">" in (67), then the three roots of Qa are found at {xi,nyg^,;, with 

Oyp-i^M < xi < 0.jg,n- Finally, if there is "=" in (67), then the roots of Qa are found at 

— kAjg,n/ l,n}- 

C. If jo = k — 1 , then the three roots of Qa are found at {Clk- 2 ,n, k^k-i,n, ^ 3 } with < 

X3 < P^k,n- 

These, combined with the fact that the leading coefficient of Qa is negative, mean that all inequal¬ 
ities in (45) are satisfied. Thus, in view of Lemma 2, < ry^. 
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Now fix an arbitrary index 0 < jo < k — 3. The inequality < rjt+yg is trivial if rk+jo = +oo- 
Otherwise, if r^+jg < +oo in (25), then, by definition, (24) and = 0 hold. So due to 

Lemma 5, (50) holds with r := and j '.— jo. By defining 

^n(^) •— (P^k,n T ^{^k,n ^) T 1^) 

and with a and b given by (50), we see that (45a) and (45c) are true because of — Ab < 0. On 
the other hand, the Ajg[r) ■ {a,b,l)^ = 0 relation in (50) expresses the fact that the inequalities 
corresponding to j G {jo, jo + 1/^ — 1} are binding in (45b). So, just as in Case II in Step 2, we 
see that + rq,, is positive on (— oo, Qy^ ^) U {Cljg+i^n, ^k-i,n), and negative on (Qy^ „, U 

(Ojt-i^n, +oo). Thus all inequalities in (45b) hold. This means that qn satisfies conditions (45) with 
this r value, so, in view of Lemma 2, C„ < The proof of the theorem is complete. □ 

8.4 The proof of Theorem 5 

Proof. We follow the ideas of the proof given in [NK16] for the fixed-step-size case. 

Suppose to the contrary that p, k and my „ (1 <;< k) satisfy the conditions of the theorem for 
some n > k and Ki, K 2 > 1 , but for all formulae with k steps and order of accuracy p, we have 
C„{co,3,l3) = 0. Then by Lemma 2, there exists a non-zero real polynomial qn that satisfies the 
conditions (45a)-(45c) with r = 0, moreover > 0 on [0, and degq„ = p (Properties 1 and 
6 of Lemma 2). 

First we define A := max^^iQqn(x) and b max^g[o,n,,„] Wn{x)\ and we introduce the 
polynomial P{x) := q„ ■ Clk^„) — j. Then the Markov brothers' inequality (see, e.g., [NK16]) 
for the first derivative implies that 

max \P'{x)\<p^- max |P(^)|, 

that is, b^^ < p^4- ^ other hand, summing the lower estimates in (32) we get ^ < Cl^^nf 
implying | ■ ^ < ^Cik,n- Thus 

bk/Ki < p^A. (68) 

Now, because of qn{Ci.k,n) = 0/ the Newton-Leibniz formula and elementary estimates yield 
that ^ 

A< [ max{0,—q'„{t))dt. (69) 

Jo 

Here we notice that the polynomial q'n is of degree at most p — 1, and q'njClj^n) > 0 for all 
0 < / < k — 1 , so the set 

{/ G Z n [0,k — 1] : G [Qy^„, ny+i^„] with q'„{x) < 0} 
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has at most p/2 elements. Therefore—^by decomposing the interval [0, as the union of the 
appropriate subintervals of length C 0 j^„, and applying the upper estimate in (32) and the estimate 
—q'n<b at most p/2 times—^we get that max( 0 , —q'^{t))dt < jK 2 b, and hence 

A < bpK2/2. (70) 

Inequalities ( 68 ) and (70) imply that k < p^K\K2/2, which contradicts the assumption of 
Theorem 5. Hence there is a formula with k steps, order of accuracy p, and C„ {co, 3, j 6 ) > 0. □ 


9 The proofs of the theorems in Section 4 

In Section 9.1 we first prove a theorem about the convergence of some rational recursions. This 
Theorem 11 will then be used in Sections 9.2 and 9.3 to prove Theorems 6 and 7, respectively. 


9.1 Global attractivity in a class of higher-order rational recursions 

Theorem 11. Let us fix an integer k > 3 and a real number A > 0. Suppose that for 1 < j < k — 1 the 
initial values Ty > 0 are given such that jj-ll Tj > 0. For any n >kwe define 


Tn ■- 


yk-l 

‘■n-] 

A + Lplrn-i 


(71) 


Then 


lim Tn 

n—>-+00 


0 , 

k-l-A 

k-1 


if k — l<A, 
if 0<A<k-l. 


To prove Theorem 11, we will apply the following lemma. 

Lemma 6 (Theorem A.0.1 in [KL02]). Suppose that a < b are given real numbers, k > 3 is a fixed 
integer, and the numbers Xj are chosen such that Xj E [a, b] for 1 < j < k — 1. Assume further that 

1 . f: [a,bf~^ —> [a,b] is continuous, 

2 . f is non-decreasing in each of its arguments, 

3. there is a unique x E [a, b] such that f{x,x,... ,x) = x, 

4. and the sequence is defined for n>k as 

■ f (^n— 1 / • • • / (/c— 1 )) • 


Then lim„^_|_oo Xn = x. 
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The straightforward proof of Lemma 6 is found in [KL02], see their Theorem 1.4.8 (for k = 3), 
Theorem A.0.1 (for k = A) or Theorem A.0.9 (for general k )—the idea of the proof is the same in 
the easiest case when / is non-decreasing in each of its arguments. Notice that in [KL02, Theorem 
1.4.8] one should have "The equation f{x,x) = x has a unique solution in [fl, fcj" instead of "... 
a unique positive solution". Now we give the proof of Theorem 11. 


The proof of Theorem 11. For some a < b {to be specified soon) we set 

Tr-l Z ■ 

/(zi.Z2.....z,_i) := 

^ + Lj=i Zj 

with zy e [a, b] (j = 1,2,... ,k — 1). Then 

A 

{djf){zi,Z 2 ,...,Zk-i) = - -^ > 0, 


(72) 


(73) 


hence / is non-decreasing in each of its arguments (and trivially continuous). Notice that due to 
(71) we have t„ G (0,1) for any n > fc, so by shifting the indices we can assume that Tj e (0,1) 
for 1 < / < A: — 1, and t„ — /(t„_i, t„- 2 , ■ ■ ■, for n >k. We distinguish two cases. 

1. The case k — 1 < A. Then / maps [a,bf^^ to [a,b] with a := 0 and b := 1. Now for any 
X E [0,1] we have 


X-f{x,x,...,x) 


x[A — {k — 1) + {k — T)x\ 
A + {k— T)x 


(74) 


precisely if x = 0, so Lemma 6 yields lim„^+oo = 0. 

2. The case 0 < A < k — 1. We set r* := mini<y<;c_i Ty. Then r* > 0 and Ty G [t*, 1] for 
1 < ; < A - 1. 

(a) The case t* > We choose a := < 1 =: b, and notice that f{a,a,... ,a) = a 

and f{b, b,... ,b) < b. By also using the non-decreasing property of / in each of its 
arguments we obtain that / maps [a, bf^^ to [a, b]. Now for any x G [a, b] the equality 
(74) holds if and only ifx = a, so Lemma 6 yields lim„^_|_oo t„ = ■ 

(b) The case 0 < r* < This time we choose a := r* < 1 =: b. Since now 

a < f{a,a,...,a) <?=> t * < f{T*,T*,...,T*) <?=> t * < ^ ' 

we have just as before that / maps [a,bf^^ to [a,b]. For any x G [a,b] we have (74) 
precisely if T = therefore we can use Lemma 6 again to get lim„^+oo t„ = 

k-l-A 
k-1 ■ 


36 












□ 


Example 1. The sequence {Tn)n>i defined by (71) can have long, non-monotonic starting slices. Consider, 
for example, the case fc = 4 with 


_ 1 _ 95638788642 

Tl , T2 — —, T3 ;^00000000000' 


and 


Tn 


Tn-l + Tn-2 + 


for n > 4. 


1 + Tn-l + + Tn-3 

T/zen 1 /ze consecutive monotone non-increasing subsequences ofT„for 1 < n < 1000 has lengths 

{2,3,2,1,2,1,2,1,2,1,2,3,3,3,3,2,1,2,1,2,1,2,1,2,3,3,3,3,2,1,2,1,2,1,2,3,3,3,3,3, 917). 


9.2 The proof of Theorem 6 

Proof. We prove the theorem for k = 3 first. We define two sequences 

+ h^-.= hr,hf-.= h2, (75) 


K-l + ^n-l + h 


and their scaled counterparts t“ := h:^ /y^, := hf /(n > 1). Then t“ and satisfy 

\-2 + \-l 


\-2 + T 2 -I + 1 


T+ = > 0, T 2 ± > 0. 


\-2 + '^n-1 + 1 


By applying Theorem 11 with k = 3 and A = 1 we see that t„ —> 1/2 hence t fz /2 as 
n —>• + 00 . Similarly, we get hf —>• y^ 12. We now define 


( 0 ,+oo)^ 9 {a,x,y) i-> f{a,x,y) := a 


x + y 


fl + X + y 

(cf. (72)). It is elementary to see that for any {a, x,y) G (0, + 00 )^ we have 
~ (x -\- 1 /)^ 

dif{a,x,y) = ,2 > ^ 2 f{a,x,y) = 83 /( 0 , x,y) = 


{a + x + yf 


[a + x + yf 


> 0 


(76) 


(77) 


(cf. (73), and notice that the function a 1 —>• f{a,x,y)/a, for example, would be monotone decreas¬ 
ing). Clearly, for zz = 1,2 we have 

K <hn< K, (78) 

so we can suppose that (78) has already been proved up to some n >2. Then by repeatedly using 
the inequality y^ < yn < y^ (implied by the assumption (33)), (77) and (78), we obtain 


^n+l — 


K-1 + K 

^n-l y^ 


- ^ hn-1 3-h _ /z„_i +hn _ 

h <T - , M < 


hn-1 -\-hn -\- y 


hn-1 -\-hn-\- y~ 
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Figure 4: The first 200 terms of the sequences hn and with }i„ := 1 + sin(n)/2, := 1/2 

and := 3/2, see (35) and (75). 



hn—\ ~\~ hfi 


hn — 1 T 




hn-1 + h„ + 


hn-1 +hn + }ln+l 

hn-1 + 


F«+l — ^n+l ^ 

K-i + K 


hn-l -\-hn + 


■r< 


K-i +hi+^+ 




= 

— “n+l- 


This shows the validity of (78) for all n > 1 by induction. By taking liminf and limsup in (78), 
Theorem 6 for k = 3 is proved. 

The proof of Theorem 6 in the general case requires only formal modifications of the argument 
given above: Theorem 11 with a general k > 3 and with A = 1 implies that for the corresponding 
sequences we have —)■ as n —)■ +oo, and the corresponding function /: (0,+oo)*^ —> 

(0, +oo) is increasing in each of its arguments. □ 


Figure 4 gives a graphical illustration of Theorem 6 for k = 3, using a hypothetical sequence 
of values for 


9.3 The proof of Theorem 7 

Proof. Step 1. Initially we suppose that some values of p > 0 and 0 < (?fe < 1 have already been 
chosen; we will make an actual choice for them in Step 5 so that the inductive argument given 
below becomes valid. First let us set n = k. 

Step 2. We know from (40) that 


hm < Q ■ for m = 1,2 ,...,n — 1. 


(79) 
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We will prove that 


Jc-l 


y"! hn—j < y/s jiy, 


and 


K < Q-hpEiUn). 

Step 3. Applying (79) and (34) repeatedly we have 

k-l k-1 Jc-l 1. \ Jc-l _ 

Yj ^n-j < Q-Y2 ^FE (««-;•) < Q'YI ^ ^Fe(m?!-i) ' 








jti (^fe)^' 


;-i' 


(80) 

(81) 

(82) 


On the other hand, the definition of in (14), the repeated application of (34), and 0 < (Ife < 1 
imply that 


\/8 > \/8^^in ^ (^((Ife)-^ ■ %E(wn-i)^ = %E(Mn-i) ' '/s (qfe)^ 


(83) 


By comparing the right-hand sides of (82) and (83) after a division by hEE(u„-i) > }i >0, we 
get that if q and Pfe are chosen such that 


J:-l 


Yh 7—TPT — (^fe)^ (? > 0, 0 < (Ife < 1/ 

(Ofe)^ 


(84) 


then (80) holds for this particular n value. 

Step 4. In order to show (81), we first notice that the function 

(0, +oo)2 9 (z, a) ^ f{z, a) := ■ a 

z + 2fl 

is increasing in each of its arguments (cf. (76)-(77)). This monotonicity property, the definition of 
hn in (39), (82), and the inequality < hEE{un-i) yield that 


hr, = 


yk-i ^ . 

4^/=i "n-i 




■ hn ^ 


%e(m«-i) ■ LLi 


■'i'—1 {qfeV 


hEE{Un-l) ■ Lyl (gJ);-i ) +2/lFE(Wn-l) 


• hEE{Un-l) 


Y^k—1 

2^/=1 


(cfeV ^ h (ii i 
7 r^l -^ • ^FE(Wn-l)- 

On the other hand, from (34) we see that q ■ //FE(wn) > Q ■ (Iee • ^EE(wn-i). Therefore if q and (Ife 
are chosen such that 

V^/c—1 Q 

T~ rrr 

< (? • (?FE/ (85) 


'1-1 (Cfe)^ ^ 


2 + 


'1=1 (m)’ ^ 
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Figure 5: The solution set of the inequalities (84)-(85). In the k = A case, the blue region has 
vertices at (1/3,1), (\/8/3,1), ~ (0.557,0.878). In the k = 5 case, the orange region has vertices 
at (1/2,1), (1/ \/2,l), « (0.569,0.9615). The points corresponding to (41) have been circled. 


then (81) holds for the actual n value. 

Step 5. So (80)-(81) will be proved as soon as we have found some and ppE satisfying (84) 
and (85). Figure 5 depicts the solution set of this system of inequalities (84)-(85) in the variables 
(q, (?fe)- The choice made in (41) is a simple rational pair; clearly, one could for example relax the 
assumption on p (by choosing it larger), but then condition (34) would in general become more 
stringent. 

Step 6 . According to the discussion preceding Theorem 7, on the one hand we have > 

2 hence C„ > 0. On the other hand, (80) guarantees h„ = Cnjin > 0 and this is the maximum 
value of hn preserving the SSP property. 

Step 7. Now we repeat Steps 2-6 inductively for each n >k + \ (Step 5 is no longer needed 
since ((?, ()fe) have already been given some particular values). The range of m in the induction 
hypothesis (79) is extended step-by-step by (81). 

Step 8 . Finally, to prove (42), we make use of the fact that 

, _ EK Zj 

{0,+Oof 3 (fl,Zi,Z2,. . .,Z)c-l) /(fl,Zi,Z2,. . .,Z;t-l) := « ' / , [ \ - 

{Z)z]zj)+2a 

is increasing in each of its arguments (cf. (76)-(77)), and repeat the steps presented in Section 
9.2: we define the corresponding auxiliary sequences and and apply Theorem 11 with 
k G {4,5} and A = 2 to show that —)■ as n —> -|-oo. □ 
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